Stochastic thermostats: comparison of local and global schemes 
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We show that a recently introduced stochastic thermostat [J. Chem. Phys. 126, 014101 (2007)] 
can be considered as a global version of the Langevin thermostat. We compare the global scheme 
and the local one (Langevin) from a formal point of view and through practical calculations on a 
model Lennard- Jones liquid. At variance with the local scheme, the global thermostat preserves the 
dynamical properties for a wide range of coupling parameters, and allows for a faster sampling of 
the phase-space. 



The most common approaches to isothermal molecular 
dynamics are perhaps those based on the introduction of 
an extended Lagrangian. The root of all these schemes 
is the Nose algorithm often used in the Hoover for- 
mulation 0. This scheme can be rigorously shown to 
provide the correct Boltzmann distribution and has a 
conserved quantity, which can be used to check the inte- 
gration timestep. A major drawback of the Nose-Hoover 
method is that it is not ergodic in some difficult cases, 
such as harmonic systems. Several different extensions 
of the Nose-Hoover method have been introduced, the 
most notable one being the so-called Nose-Hoover chains 
Q, which addresses the ergodicity issue at the price of 
an increased complexity in the algorithm. Although the 
Nose-Hoover scheme was originally written as a global 
thermostat, i.e. coupled only to the total kinetic energy 
of the system, it is sometimes implemented in a local 
manner (also called massive Nose- Hoover) , i.e. using an 
independent thermostat on each degree of freedom |J| . 

Another common choice is the weak-coupling method, 
introduced by Berendsen et al [f|. This scheme is a con- 
tinuous version of the velocity-rescaling scheme, thus it is 
a global thermostat. It is deterministic, stable and easy 
to implement, but it does not produce configurations in 
the canonical ensemble. 

An alternative approach to canonical sampling is to use 
stochastic molecular dynamics. The most common form 
is Langevin dynamics [6( . The Langevin thermostat is lo- 
cal, and its major feature is that ergodicity can be proven 
also in pathological cases. However, since the friction and 
noise terms alter significantly the Hamiltonian dynamics, 
it cannot be used to compute dynamical properties, un- 
less an extremely small friction is used. Moreover, the 
effect of the friction and noise terms on the sampling 
efficiency is non trivial. Even in applications where dy- 
namical properties are not relevant it can be difficult to 
properly tunc the friction in order to achieve an efficient 
sampling. 

In a recent paper 0] we proposed a stochastic veloc- 
ity rescaling which can be considered as Berendsen ther- 
mostat plus a stochastic correction leading to canonical 



sampling. We also showed that, in spite of its stochastic 
nature, one can define a conserved quantity. This scheme 
does not suffer of ergodicity problems in solids Q, has 
been used in practical applications for equilibration pur- 
poses [|[ or to perform ensemble averages [H, [To[ and can 
be combined with variable-cell dynamics to perform sim- 
ulations in the isothermal- isobaric ensemble In the 
present paper we present an alternative derivation of the 
same scheme, where stochastic velocity rescaling is ob- 
tained starting from Langevin dynamics and minimizing 
the disturbance of the thermostat on the Hamiltonian 
trajectory, nevertheless retaining the same thermaliza- 
tion speed of Langevin dynamics. This idea was also 
used by Berendsen et al to derive their algorithm [f|. 
Moreover, we show how stochastic velocity rescaling can 
be considered as a global version of Langevin dynamics. 
Thus the relationship between the two schemes is similar 
to that between standard Nose-Hoover and massive Nose- 
Hoover. Finally, we compare in practical situations the 
efficiency of the local (Langevin) and global (rescaling) 
versions, and show that the disruption of Hamiltonian 
dynamics observed using Langevin thermostat is not due 
to the the stochastic nature of the algorithm but to the 
use of a local thermostat. 



I. CONTINUOUS EQUATIONS OF MOTION 

We consider a system described by coordinates and 
momenta p,-, where i runs over the Nf degrees of freedom, 
and with q and p wc indicate the set of coordinates qi and 
Pi. We associate a mass m.j to each degree of freedom, 
and we define a Hamiltonian H(p,q) = K{p) + U(q), 
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where U(q) is the potential energy, and K(p) = J2i 
is the kinetic energy. We want to sample the canonical 
distribution P(p, q)dpdq oc e^C^H^M) , where (3 is the 
inverse temperature, by means of equations of motion in 
the form 



, / \ dU , , , , 
dPi(t) = -—dt + gi(t)dt 

dqi (t) = *®dt. 
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Equations ([T]) are Hamilton equations plus a correction 
force gi (t) which artificially modifies the dynamics of the 
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FIG. 1: Schematic representation of the momentum compo- 
nents at time t and at time t + At. g(t) is a generic force 
applied to the system. g(t) is a force which leads to the same 
change in the kinetic energy as g(t), but minimizes the dis- 
turbance. 



system. Since the total energy H is conserved in Hamilto- 
nian dynamics, only gi(t) is responsible for its variations 
and leads to the system thermalization. 

In standard Langevin dynamics, the correction force is 



gi{t)dt = -jp^dt 



(2) 



where 7 is the friction coefficient, and dWi(t) is a vec- 
tor of Nf independent Wiener noises, normalized as 

( dW dt t} dW dt' ] ) =<*(*- Thc thermalization speed 

can be quantified calculating the time derivative of the 
total energy from Eqs. ((T|) and @. Using thc Itoh chain 
rule [121 we obtain 



dH ( t )=y \-^dt+ 



'27 PM 
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dW.it) 



+ ^*. (3) 

This expression can be further simplified defining the av- 
erage kinetic energy K = iV/(2/?) _1 , a relaxation time 
r = (27) -1 , and exploiting the fact that the noise terms 
on different degrees of freedom are independent of each 
other: 



dH(t) 




dW(t). (4) 



but minimizes thc disturbance on the trajectory. This 
procedure is exactly the same used by Berendsen et al 
Q , the only difference being that Eq. ^ on the total en- 
ergy now contains a stochastic term. We first notice that, 
since the force only acts on the momenta, fixing a value 
for H is equivalent to fixing a value for K . Following 
Ref. Q we quantify the disturbance as J ^ li m~ 1 igiit)dt) 2 . 
As it is seen in Fig. [TJ the minimal disturbance for a fixed 
kinetic energy increment is obtained with a force <?i(i) 
which is proportional to Piit). Thus ffj(t) = A(i)pj(t), 
where A(t) is chosen so as to enforce a given variation of 
the total energy. Since A(i) includes a stochastic part, 
and since the variation of the total energy depends on 
the momenta only through the kinetic energy K, this 
last relation can be written as 

gi(t)dt = Pi (t) [A(K(t))dt + B(K(t))dW(t)] , (5) 

where AiK) and -B(-ftT) are arbitrary functions of the ki- 
netic energy. Thc change in the total energy is then 

dH = 2A(K)Kdt + 2BiK)KdW + B 2 iK)Kdt. (6) 

Expressions for AiK) and B[K) can be obtained setting 
Eq. (0| equal to Eq. ©, resulting in the correction force 



g%dt 
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N f Kr 



P.dW. (7) 



This equation is stochastic, with the same noise term 
used on all the particles. It is also very similar to the 
expression of thc force in the Berendsen algorithm. 

The combination of Eqs. |T]) and ([7]) results in a contin- 
uous, stochastic dynamics which can be shown to sample 
exactly the canonical ensemble. The effect of a gt par- 
allel to pi is the same of a rescaling procedure and the 
enforced increment of the total energy in Eq. (g]) is the 
same as in Ref. Q. Thus, Eqs. ([I]) and ([7]) represents the 
continuous version of thc velocity rescaling described in 
Ref. Q. 

Notably, if Nf = 1, Eq. ([7]) becomes equivalent to 
Eq. ^j. Thus when the thermostat is applied to a sin- 
gle degree of freedom, it is completely equivalent to a 
Langevin thermostat. One can perform Langevin molec- 
ular dynamics by applying a thermostat per degree of 
freedom, or stochastic rescaling by applying a single ther- 
mostat to the total kinetic energy. Intermediate schemes 
can be designed, where a thermostat is applied on each 
molecule or group of atoms. 



It is worth noting that while in Eqs. ([TJ) and ([3]) there 
are Nf independent noise terms, in Eq. ((¥]) a single noise 
term is present. 

We now want to design a new correction force giit) 
which gives the same variation of the total energy as 
Langevin dynamics, thus the same thermalization speed, 



II. FINITE TIMESTEP ALGORITHM 

In the practical implementation, time is incremented 
in discrete steps, and the Trotter decomposition scheme 
[l3l [HI can be used to separate the integration of Hamil- 
ton equations and the update of the momenta due to g. 
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The former is then integrated using standard velocity- 
Verlet, while for the latter we need to integrate Eq. {7J. 
A possible approach is to consider the propagation of 
kinetic energy when the momenta evolve according to 
Eq. {?]), as it is done in the appendix of Ref. @- The 
analytical solution of Eq. ([7|) for a finite time At is 



Pi (t + At) =a(t) Pi (t), 



(8a) 



where 



a 2 {t) = c + 




N f K{t) 



(8b) 



Here c = e~ 2 ^ At = e~ At / r , R(t) is a Gaussian number 
with unitary variance and S^f-i is the sum of Nf — 1 
independent, squared, Gaussian numbers. Equation ([8|) 
has been obtained enforcing the evolution of the kinetic 
energy, thus, strictly speaking, it does not fix the sign 
of a. A more rigorous analysis shows that the sign of a 
should be chosen as 



sign(a(t)) = sign R(t) 



l cN f K(t) \ 
(l-c)R)' 



(9) 



to keep into account the finite probability to observe a flip 
of the momenta pi when the force in Eq. ((7]) is applied. 
The Gaussian number in Eq. ^ needs to be the same 
that is used in Eq. (|8|). The probability to observe the 
flip is extremely small if Nf is large and c « 1, which 
is the usual case when the thermostat is used as global 
and r > At. This is the case in Ref. Q, where we set 
a > 0. On the other hand, when the thermostat acts 
on a few degrees of freedom, the sign of a needs to be 
calculated by means of Eq. ©. This is always the case 
for Langevin dynamics. With simple manipulation, it 
can be shown that for Nf = 1 the integration scheme 
given by Eqs. © and ©, combined with velocity- Ver let, 
is completely equivalent to the integration scheme for 
Langevin dynamics introduced in Ref. [l4| . 



III. EXAMPLES 

Up to now we simply established a theoretical relation- 
ship between Langevin dynamics and stochastic rescal- 
ing. The outcome is that the effect of the two algorithms 
on the total energy should be equivalent if the friction 
in the Langevin dynamics and the relaxation time in 
the stochastic scale are related by r = (27) . How- 
ever, the stochastic rescaling is expected to give better 
dynamical properties, since only the component of the 
force which changes the total energy is retained. We 
test here this affirmations on a simple test-case, namely 
a Lennard-Jones fluid with density p = .8442 and tem- 
perature /3 _1 = 0.722, which is close to the triple point. 
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FIG. 2: Diffusion coefficient as obtained from thermostated 
simulations as a function of the thermostat relaxation time 
r, for the local and global thermostats as indicated. The sta- 
tistical error is smaller than the symbol size. As a reference, 
the diffusion coefficient of a free particle subject to the same 
Langevin equation is plotted (dashed line) . All the quantities 
are in Lennard-Jones reduced units. 



Throughout this section we use reduced Lennard-Jones 
units for temperature, distance and time. We simulate 
a box containing 108 particles, with periodic boundary 
conditions, and we cut the interaction at distance 2.5. 
We set the timestep to At = 0.005, which leads to a 
reasonable conservation of the effective energy 0, fh4| . 
We then perform runs of 10 7 steps, with both the local 
scheme (Langevin dynamics) and the global one (stochas- 
tic rescaling), using a broad range of values of the ther- 
mostat relaxation time r. 

To quantify the disturbance on Hamiltonian dynamics, 
in Fig. [2] we show the diffusion coefficient as a function 
of t, as obtained from the Einstein relations. When the 
local thermostat is used with a short relaxation time, 
the diffusion is strongly quenched. This happens when 
the typical collision time with the external bath 7 _1 is 
shorter than the typical collision time between the par- 
ticles, so that the former becomes the real bottleneck for 
the diffusion process. In the limit of short r, the equa- 
tions of motion tend to a high-friction Langevin dynam- 
ics. In this case, the observed diffusion coefficient D is 
proportional to (flmj)' 1 , which is the value of the dif- 
fusion coefficient for a free particle subject to the same 
Langevin dynamics. The prefactor is related to the dif- 
ficulty in crossing the barriers between different liquid 
configurations. On the contrary, with the global thermo- 
stat D is almost independent on r, indicating that the 
disturbance on the dynamics is very small and that good 
estimates of D in the canonical ensemble can be obtained 
also with a thermostated simulation. 

To quantify the equilibration speed we calculate the 
autocorrelation time of a few global observables. The 
efficiency of a sampling algorithm is optimal when the 
autocorrelation time tx of the desired quantity X is min- 
imal. If X is the total energy, tx also indicates how fast a 
simulation started from an unlikely configuration is equi- 
librated. We integrate the autocorrelation function using 
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FIG. 3: Autocorrelation time of the kinetic (tk), potential 
(tu), and total (th) energy, as a function of the thermostat 
relaxation time r, for the local and global thermostats as 
indicated. The statistical error is smaller than the symbol 
size. All the quantities are in Lennard- Jones reduced units. 



choose T = 50, and we expect a relative accuracy on tx 
on the order of 5%. Since T ^> t Xi t x is also a good 
approximation for the autocorrelation time. 

In Figure [3] we plot the autocorrelation time of the ki- 
netic energy tk, of the potential energy tjj and of the 
total energy th, as a function of the thermostat relax- 
ation time t. The autocorrelation time of the kinetic 
energy is completely dictated by the thermostat relax- 
ation time, and independent on the choice of a local or 
a global scheme. On the contrary, the autocorrelation 
time of the total energy and of the potential energies are 
proportional to r only in the limit of large r. In the local 
scheme, when r is smaller than 0.2 the disturbance of the 
trajectory becomes so large that the phase-space explo- 
ration turns out to be slower. Comparing Figs. [2] and [3j 
it is seen that the optimal value for r is the smallest one 
that still does not affect the diffusion coefficient. When 
the global scheme is adopted, even small values of r can 
be safely used, resulting in a faster decorrelation of the 
total energy. 

IV. CONCLUSIONS 



a windowing function 



T X = 



(It 



(5X(0)5X{t)) 



(6X(0)SX(0)) 



1 - 



T 



(10) 



where SX = X — (X) and T is a large value. The 
windowing function in parcthcsis helps the convergence 
because it weights less the points with larger statisti- 
cal error. Moreover, Equation (jTUJ) is exactly equal to 
Te 2 {T)/(2(6X 2 )), where e 2 (T) is the mean square er- 
ror from a time average of length T. The relative accu- 
racy in evaluation of tx is approximately y // 2T/T', where 
T" = 5 x 10 4 is the total run length. In the following we 



In conclusion, we have presented an alternative deriva- 
tion of the global thermostat introduced in Ref. Q • This 
derivation allows to write continuous equations of motion 
and shows the analogy between this scheme and the stan- 
dard Langevin thermostat. Namely, the new scheme can 
be considered as a global version of the Langevin ther- 
mostat, that minimizes the disturbance of the original 
Hamiltonian dynamics. Finally, we have discussed these 
properties on a simple test case, showing that the global 
scheme preserves the dynamical properties. Moreover, 
using as a measure the autocorrelation time of the total 
energy, we have shown that the global scheme allows for 
a faster sampling. 
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